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Abstract.  Computer  program  HEC-1 ,  a  precipitation-runoff  model  widely  used 
throughout  the  United  States,  includes  the  capability  to  estimate  automati¬ 
cally  any  of  twelve  parameters  necessary  to  model  the  precipitation-runoff 
process  and  the  channel  routing  process.  The  parameter  estimation  scheme 
employs  Newton's  method  to  minimize  a  weighted  sum  of  squares  of  differences 
between  observed  and  computed  hydrograph  values.  Applications  of  this  pa¬ 
rameter  estimation  procedure  are  presented,  and  typical  steps  of  the  procedure 
for  determining  optimal  parameter  estimates  are  outlined.  Recent  efforts  to 
Improve  the  estimation  algorithm  and  recent  use  of  the  calibration  capability 
to  update  sequentially  parameter  estimates  in  a  flood  forecasting  application 
are  discussed. 

Keywords .  Water  Resources,  Hydrology,  Rainfall -Runoff  Modelling,  Parameter 
Estimation,  Computer  Applications. 


INTRODUCTION 

Computer  program  HEC-1,  which  simulates  the 
hydrologic  response  of  urban  and  rural  water¬ 
sheds,  was  developed  as  a  tool  to  assist  the 
staff  of  the  U.S.  Army  Corps  of  Engineers  in 
meeting  their  water  management  responsibili¬ 
ties.  The  basic  concepts  embodied  in  the  pro¬ 
gram  were  conceived  in  1966  when  Beard  and 
other  members  of  The  Hydrologic  Engineering 
Center  (HEC)  developed  a  set  of  small  programs 
that  could  be  used  independently  to  solve  the 
Individual  tasks  typically  required  in  a  hydro- 
logic  study.  Included  In  the  set  of  programs 
was  one  which  employed  a  univariate  version  of 
Newton’s  technique  to  calibrate  automatical ly 
unit  hydrograph  and  loss  rate  parameters 
(HEC,  1967).  In  1967,  when  the  set  of  programs 
was  combined  into  the  single  computer  program 
entitled  HEC-1  Flood  Hydrograph  Package,  this 
technique  was  adopted  and  has  been  retained 
through  subsequent  revisions,  including  the 
latest  1973  version.  Over  400  copies  of  the 
latest  version  have  been  distributed  to  pri¬ 
vate  consulting  firms,  universities,  and 
governmental  agencies  in  both  the  United 
States  and  other  countries.  The  program  has 
been  executed  nearly  4000  times  annually  by 
Corps  personnel  using  the  Lawrence  Berkeley 


Laboratory  computer. 

The  automatic  parameter  estimation  option  of 
the  program  has  been  extensively  used  for  sev¬ 
eral  reasons.  First,  the  Hydrograph  Analysis 
Package  is  composed  of  a  set  of  simple  con¬ 
ceptual  procedures  employing  lumped  parameters 
that  are  intended  to  model  the  general  behav¬ 
ior  of  hydrologic  phenomena.  The  parameters 
of  the  conceptual  procedures  are  inferred  from 
precipitation  and  streamflow  records  rather 
than  by  direct  measurement  of  physical  water¬ 
shed  characteristics.  As  a  result  of  inac¬ 
curacies  in  the  modelling  process  and  in  the 
measurement  of  input  data,  even  experienced 
users  have  difficulty  in  determining  precise 
parameter  values.  Second,  many  users  are  in¬ 
experienced  with  the  model,  and  the  automatic 
calibration  feature  eliminates  the  frustrat¬ 
ing  trial  and  error  approach  to  estimating 
acceptable  parameter  values. 

This  paper  describes  the  major  technical  com¬ 
ponents  of  the  program  and  details  the  basic 
structure  of  its  calibration  technique.  The 
strategy  for  employing  the  automatic  calibra¬ 
tion  features  of  the  program  In  regional  stud¬ 
ies  Is  discussed,  and  several  typical  applica¬ 
tions  are  described.  A  flood  forecasting 
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application  is  suonarized  that  involves 
sequentially  updating  the  model's  parameters 
as  successive  forecasts  are  calculated  during 
a  storm  event.  In  addition,  application  of 
an  alternative  optimization  scheme,  involving 
the  random  search  method  of  Box  (1965),  is 
presented. 


DESCRIPTION  Of  THE  HEC-1  FLOOD 
HYDROGRAPH  PACKAGE 

The  HEC-1  computer  program  consists  of  three 
major  hydrologic  components  that  determine 
the  average  precipitation  and  snovmielt 
and  the  amount  of  effective  precipitation 
contributing  to  direct  runoff  from  a  subbasin, 
compute  the  subbasin  runoff  hydrograph  from 
the  effective  precipitation,  and  route  and 
combine  the  subbasin  runoff  hydrographs.  All 
components  In  HEC-1  employ  lumped  parameters 
for  each  subbasin  or  routine  reach.  This  means 
that  the  model's  input,  parameters,  and  output 
are  considered  to  be  average  values  over  the 
entire  subbasin  of  Interest.  When  subbasin 
averages  are  not  appropriate,  smaller  subbasins 
can  be  defined  to  obtain  a  better  spatial 
definition. 

The  first  computational  step  in  the  program 
determines  the  average  subbasin  precipitation 
from  either  historical  gaged  data  or  hypo¬ 
thetical  storms.  This  is  followed  by  an 
accounting  of  interception  and  accumulation 
of  soil  moisture  by  a  loss  rate  function  to 
compute  the  subarea  rainfall  or  snowmelt 
excess.  The  moisture  excess  is  then  distri¬ 
buted  in  time  by  a  unit  hydrograph  function, 
added  to  a  base  flow  function,  and  recessed 
by  a  logarithmic  decay  function  once  a  spec¬ 
ified  recession  flow  rate  is  reached  on  the 
falling  limb  of  the  hydrograph.  .  This  yields 
the  total  runoff  hydrograph  at  the  subbasin's 
outlet.  Next,  the  runoff  hydrograph  from  the 
subbasin  is  conveyed  downstream  using  a  stream- 
flow  routing  function.  Tributary  streamflows 
are  added  at  the  confluences  as  the  simulation 
proceeds  downstream.  Snowfall  and  snowmelt 
are  simulated  in  each  subarea  according  to 
temperatures  In  various  elevation  bands. 

Either  the  degree-day  or  energy  budget  method 
may  be  used  to  compute  the  snovmelt.  Complex 
stream  systems  can  be  simulated  if  the  unit 
graph,  loss  rate  and  routing  parameters  are 
specified  along  with  either  observed  precipi¬ 
tation  or  a  specified  precipitation  depth-area 
relationship.  In  addition,  expected  annual 
damages  can  be  calculated  at  any  point  in  the 
stream  system  where  flow- frequency  and  flow- 
damage  relationships  are  provided.  The  model’s 
automatic  calibration  features  can  be  used  to 
select  unit  hydrograph  and  loss  rate  parameters 
in  a  single  subbasin  or  choose  the  routing 
parameters  In  an  Individual  river  reach.  In 
both  cases  the  optimization  scheme  Is  based 
on  comparisons  with  the  observed  and  simulated 
hydrographs. 


Loss  Rate  Functions 

Precipitation  losses  to  interception,  depres¬ 
sion  storage  and  Infiltration  may  be  simulated 
by  one  of  three  loss  rate  functions:  initial 
loss  followed  by  a  constant  loss  rate,  the  SCS 
curve  number  technioue  (Soil  Conservation 
Service,  1972),  or  the  HEC  exponential  loss 
rate  function.  The  latter  computes  precip¬ 
itation  losses  as  a  function  of  antecedent 
soil  moisture,  precipi tation  intensity,  and 
an  infiltration  rate  that  is  a  non-linear 
function  of  accumulated  losses,  as  shown  in 
Fig.  1.  For  a  snow-free  basin  the  equations 
for  the  HEC-1  exponential  loss  rate  are  as 
follows: 

AlOSS  *  (AK  +  DLTK)  (RAIN)ERAIN  (1) 

AK  *  (STRKR)/(RTIOL)0,1  CUML  (2) 

DLTK  *  0.2  DLTKR  [MCUML/OLTKR)]2 

for  (CUML/OLTKR)  <1; 

otherwise  DLTK  3  0.  (3) 

where  ALOSS  is  the  loss  rate  in  nin/hr;  AK  is 
the  basic  loss  coefficient,  DLTK  is  the  incre¬ 
mental  loss  coefficient;  RAIN  is  the  rainfall 
intensity  in  mm/hr;  ERAIN  is  an  exponent  that 
reflects  the  influence  of  precipitation  inten¬ 
sity  on  the  basin  average  loss  characteris¬ 
tics;  STRKR  is  the  starting  value  of  the  basic 
loss  index  on  the  exponential  recession  curve 
in  m/hr;  RTIOL  is  the  ratio  of  the  loss  rate 
coefficient  AK  to  that  after  254  mm  more  of 
accumulated  loss  occurs;  CUML  is  the  accumu¬ 
lated  loss  in  mm;  and  01TKR  is  an  incremental 
loss  index. 

The  program  contains  a  separate  set  of  loss 
rate  equations  that  are  employed  when  the 
snowmelt  capabilities  of  the  program  are 
desired: 

ALOSS  =  AK  (RAIN  +  SNWMT)ERAIN  (4) 

AK  *  (STRKSJ/CRTIOK)0-1  CUML  (5) 


where  SNWHT  is  the  snowmelt  in  mm/hr,  STRKS 
is  the  basic  loss  coefficient  for  snowmelt 
in  ron/hr,  and  RTI0K  is  similar  to  RTIOR  for 
snowmelt  conditions.  Equations  4  and  5  are 
used  in  lieu  of  equations  1  and  2  whenever 
snowmelt  occurs.  The  amount  of  snowmelt  is 
calculated  separately  in  each  elevation  zone 
based  on  the  air  temperature  which  is  cal¬ 
culated  from  a  base  temperature  at  the  lowest 
elevation  and  a  user  supplied  adiabatic  lapse 
rate.  The  degree-day  method  for  computing 
sncnwnelt  is: 

SNWMT  *  COEF  (TMPR  -  FRZTP)  (6) 

where  TMPR  is  the  air  temperature  in  °r.  lapsed 
to  the  midpoint  of  the  elevation  zon*»,  fR7^P 
is  the  temperature  in  9C  at  which  snow  melts, 
and  COEF  is  the  melt  coefficient  in  m  per 
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degree-day  (°C).  Energy-budget  equations 
for  melt  during  rain  or  melt  during  rainfree 
periods  can  also  be  employed.  The  losses 
are  subtracted  from  rainfall  and  snowmelt  in 
each  zone,  and  the  excesses  are  sunned  to 
yield  the  excess  precipitation  from  the 
subbasin. 

In  the  loss  rate  and  snowmelt  equation,  the 
following  parameters  must  be  determined  by 
calibration;  STRKR,  RTIOL,  DLTKR,  ERAIN,  COEF, 
STRKS,  RTIOK,  and  FRZTP. 

Unit  Hydrograph  Functions 

The  precipitation  excess -to- runoff  transfor¬ 
mation  Is  accomplished  by  the  use  of  the  unit 
hydrograph.  Sheman  (1932)  defined  the  unit 
hydrograph  as  follows: 

If  a  given  one-day  rainfall  produces  a 
one- inch  depth  of  runoff  over  the 
given  drainage  area,  the  hydrograph 
showing  the  rates  at  which  the  runoff 
occurred  can  be  considered  a  unit 
graph  for  the  watershed. 

Application  of  this  technique  to  precipitation 
excess  amounts  other  than  one  inch  is  accom¬ 
plished  by  multiplying  precipitation  excess 
amounts  by  the  unit  hydrograph  ordinates 
because  the  runoff  ordinates  for  a  given  dura¬ 
tion  are  assumed  to  be  directly  proportional 
to  the  rainfall  excess.  Unit  hydrograph  ordin¬ 
ates  can  be  supplied  directly  to  the  program 
or  the  unit  hydrograph  can  be  calculated  using 
techniques  proposed  by  Clark  (1945),  Snyder 
(1938)  or  the  SCS  (1972). 

The  Clark  method  uses  two  parameters  and  a 
time-area  relationship  to  define  an  instan¬ 
taneous  unit  hydrograph.  Experience  by  the 
HEC  has  indicated  that  the  use  of  a  detailed 
time  area  relationship  is  usually  not  warran¬ 
ted  and  that  one  based  on  a  generalized  water¬ 
shed  shape  (contained  within  the  program)  is 
satisfactory  in  most  instances.  When  used  in 
this  fashion,  Clark's  technqiue  requires  only 
TC,  the  time  of  concentration,  and  R,  a 
storage  constant,  to  define  the  ordinates  of 
the  unit  hydrograph.  This  function  is  attrac¬ 
tive  because  it  avoids  the  difficulties  assoc¬ 
iated  with  the  calibration  of  many  individual 
unit  graph  ordinates.  The  general  shape  of 
the  hydrograph  is  fixed  and  problems  assoc¬ 
iated  with  negative  ordinates  and  infeasible 
fluctuations  of  the  unit  hydrograph  are  elim¬ 
inated. 

Tn  addition  to  the  Clark  parameters,  Snyder's 
coefficients  TP  and  CP,  which  define  the  peak 
of  the  unit  hydrograph,  can  be  provided  as  in¬ 
put.  The  program  internally  uses  the  Clark 
procedure  by  interatively  varying  TC  and  R 
until  the  peak  of  the  unit  hydrograph  corres¬ 
ponds  to  the  one  described  by  the  specified 
Snyder's  parameters.  The  SCS  dimensionless 
unit  hydrograph  technique,  which  uses  a  single 
lag  parameter  TLAG  to  define  the  shape  of  a 
triangular  unit  hydrograph,  can  also  be  used. 


The  unit  hydrograph  parameters  required  for 
calibration  are  TC  and  R,  or  TP  and  CP,  or 
TLAG. 

Streamflow  Routing 

Subbasin  runoff  is  routed  downstream  using  one 
of  the  following  'hydrologic*  routing  tech¬ 
niques:  Muskingum,  working  RAO,  straddle- 
stagger,  Tatum,  modified  Puls  or  multiple 
storage.  Parameters  for  the  first  four  meth¬ 
ods  can  be  automatically  calibrated.  These 
parameters  include  the  following: 

1.  Humber  of  routing  steps  to  be  used  for 
routing  by  the  Tatum  method,  Muskingum 
method,  or  modified  Puls  method  (NSTPS). 

2.  Humber  of  ordinates  to  be  averaged  in  the 
straddle-stagger  routing  (NSTDl). 

3.  Number  of  intervals  the  hydrograph  is  to 
be  lagged  in  the  straddle-stagger  routinq 
(LAG). 

4.  Coefficients  of  the  Muskingum  routing 
function  (AMSKK  and  X). 

5.  Time-of- storage  coefficient  for  the 
multiple  storage  routing  (TSK). 

The  non-linear  storage  outflow  relationship 
required  for  the  modified  Puls  and  working 
RAD  methods  can  be  obtained  using  steady- 
state  water  surface  profile  computations  or 
using  a  separate  optimization  program  such 
as  the  one  suggested  by  Slocum  and  Dandekar 
(1975). 

PARAMETER  ESTIMATION  TECHNIQUES 

If  HEC-1  were  a  perfect  model  of  watershed 
hydrology,  and  if  total  precipitation  and 
total  direct  runoff  could  be  measured 
accurately,  the  parameters  of  the  precipita¬ 
tion-runoff  transformation  functions  for  a 
particular  storm  event  could  be  determined 
directly  by  inverse  solution  of  the  transfor¬ 
mation  equations.  However,  these  conditions 
are  not  satisfied  in  reality,  and  inverse 
solution  of  the  equations  is  difficult,  so 
the  parameters  are  found  instead  by  selec¬ 
tion  of  those  values  that  yield  the  "best" 
reproduction  of  a  measured  runoff  event  with 
the  available  measured  precipitation  data  and 
the  available  model.  This  parameter  selec¬ 
tion  has  often  been  accomplished  by  a  system¬ 
atic  trial-and-error  procedure:  parameter 
values  are  selected,  the  model  is  exercised 
with  these  values,  and,  the  resulting  runoff 
hydrograph  is  compared  with  the  observed 
hydrograph.  If  the  "fit"  is  less  than  satis¬ 
factory,  different  parameter  values  are  se¬ 
lected,  the  model  is  exercised  with  these 
values,  and,  the  resulting  runoff  hydrograph 
is  compared  with  the  observed  hydrograph.  If 
the  "fit"  is  less  than  satisfactory,  different 
parameter  values  are  selected,  and  the  entire 
process  is  repeated. 


An  alternative  to  the  trial-and-error  approach 
to  parameter  selection  Is  an  automatic  cal¬ 
ibration  approach  in  which  the  necessary  tasks 
for  calibration  are  automated.  Automatic  cal¬ 
ibration  requires  selection  of  an  explicit 
Index  of  the  acceptability  of  alternative  par¬ 
ameter  estimates,  definition  of  the  range  of 
feasible  values  of  the  parameters,  and  develop¬ 
ment  of  some  techniques  for  correction  of  the 
parameter  estimates  until  the  “best*1  estimates 
are  determined.  Thus  the  parameter  estimation 
problem  can  be  classified  as  an  optimization 
problem:  there  is  an  objective  function  for 
which  an  optimal  value  is  sought,  subject  to 
certain  constraints  on  the  decision  variables 
(the  parameters).  The  HEC-1  program  includes 
the  capability  to  solve  this  optimization  pro¬ 
blem,  thereby  automatically  determining  op¬ 
timal  parameter  estimates. 

Objective  Function 

The  objective  function  of  the  parameter  esti¬ 
mation  optimization  problem  must  define  the 
difference  between  the  computed  runoff  hydro¬ 
graph  (with  any  parameter  estimates)  and  the 
recorded  runoff  hydrograph.  Presumably,  this 
difference  will  be  at  a  minimum  for  the  opti¬ 
mal  parameter  estimates.  In  HEC-1,  the  fol¬ 
lowing  objective  function  is  employed  as  an 
Index  of  the  errors: 


STDER  -7  z  [(Q0BSi  -  QCOMP,)2*  WT  J  /  N  (7) 
1*1  71 


These  constraints  are  summarized  In  Table  1. 
The  constraints  shown  here  are  limited  to 
those  imposed  explicity  In  the  program. 
Additional  constraints  may  be  appropriate  In 
certain  circumstances;  however,  these  must  be 
Imposed  externally  to  the  program  when  the 
user  must  decide  whether  to  accept,  to  modify, 
or  to  reject  a  given  parameter  set,  based  on 
engineering  judgment. 

TABLE  1  Constraints  on  HEC-1 

Uni  Oiydroqraph  and  Loss  llate 

Parameters 

(TC  +  R)  >  1. 03/(1.  -  [R/(TC  +  R)]) 

R/(TC  +  Rj  <  0 . 52/ (TC  +  R) 

R/(TC  ♦  R)  <,  1.0  -  (1.03/(TC  +  R)) 

ERAIN  <_  1.0 
RTIOL  i  1.0 
RTI0K  >  1.0 

FRZTP  7-1.11  °C 
FRZTP  <  3.33  °C 


OPTIMIZATION  TECHNIQUE 

The  constrained  optimization  scheme  employed 
in  HEC-1  is  a  univariate  search  technique 
that  uses  Newton's  method.  Application  of 
such  a  technique  permits  use  of  the  simula¬ 
tion  capabilities  of  HEC-1  in  a  traditional 
manner  and  does  not  require  development  of 
analytical  derivatives.  Steps  in  applica¬ 
tion  of  this  technique,  as  implemented  in 
HEC-1,  are  as  follows: 


where  STDER  *  the  error  index;  QOBS.  * 
observed  runoff  hydrograph  ordinateTfor  time 
period  i;  QC0MP.  =  the  runoff  hydrograph 
ordinate  for  time  period  1,  computed  by  HEC-1 
with  the  current  parameter  estimates;  N  * 
total  number  of  hydrograph  ordinates;  WT-j  * 
a  weight  for  the  hydrograph  ordinate.  The 
weight,  WT^,  is  defined  as  follows: 

WT1  *  (Q0BS1  +  QAVE)  /  (2  *  QAVE)  (8) 

where  QAVE  *  the  average  computed  discharge. 
This  weighting  function  emphasizes  accurate 
reproduction  of  peak  flows  rather  than  low 
flows  by  biasing  the  objective  function.  Any 
errors  for  discharge  ordinates  that  exceed 
the  average  discharge  will  be  weighted  more 
heavily,  and  hence  the  optimization  scheme 
should  focus  on  reduction  of  these  errors. 

Constraints 

The  range  of  feasible  values  of  the  parame¬ 
ters  1$  bounded  because  of  physical  limita¬ 
tions  on  the  values  that  the  various  unit 
hydrograph,  loss  rate,  and  snowmelt  parame¬ 
ters  may  have,  and  also  because  of  numerical 
limitations  Imposed  by  the  mathematical  func¬ 
tions  employed  to  model  watershed  behavior. 

In  addition  to  bounds  on  the  maximum  and  min¬ 
imum  values  of  certain  parameters,  the  inter¬ 
action  of  some  parameters  is  also  restricted 
because  of  physical  or  numerical  limitations. 


1.  Initial  values  are  assigned  for  all  par¬ 
ameters.  These  values  may  be  assigned 
by  the  program  user,  or  program-assigned 
default  values  may  be  used.  The  default 
parameter  values  are  shown  in  Tables 

2  and  3. 

2.  The  response  of  the  watershed  is  simula¬ 
ted  with  the  initial  parameter  estimates, 
and  the  value  of  the  objective  function 
is  computed  by  comparison  of  the  ordin¬ 
ates  of  the  computed  and  observed  runoff 
hydrographs. 

3.  In  the  order  shown  in  Tables  2  and  3, 
each  parameter  to  be  estimated  is  de¬ 
creased  by  one  percent  and  then  by  two 
percent,  the  system  response  is  eval¬ 
uated,  and  the  objective  function  cal¬ 
culated  for  each  change,  respectively. 
This  gives  three  separate  system  evalu¬ 
ations  at  equally-spaced  values  of  the 
parameter  with  all  other  parameters  held 
constant.  The  "best"  value  of  the  par¬ 
ameter  is  then  estimated  using  Newton's 
method. 

4.  Step  3  is  repeated,  using  the  "best" 
estimates  of  the  parameters. 

5.  Step  3  is  repeated  for  the  parameter  that 
most  improved  the  value  of  the  objective 
function  In  its  last  change  until  no  sin¬ 
gle  change  in  any  parameter  yields  a 
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reduction  of  the  objective  function  of 
more  than  one  percent. 

6.  One  more  complete  search  of  all  param¬ 
eters  is  made. 

7.  Step  5  Is  repeated,  and  the  final  param¬ 
eter  estimates  are  identified  as  optimal. 

For  the  Integer- valued  parameters,  the 
estimate  Is  increased  or  decreased  for  each 
in  turn  until  a  minimum  1$  determined. 

The  scheme  employed  for  estimating  the  "best" 
value  of  each  parameter  in  Step  3  is  based 
on  the  concept  that  the  optimum  of  the  objec¬ 
tive  function  occurs  at  a  root  of  the  first 
partial  derivative  of  the  function  with 
respect  to  each  of  the  parameters.  These 
derivatives  cannot  easily  be  evaluated 
analytically  because  the  objective  function 
indirectly  Includes  all  the  functions  and 
equations  contained  in  the  HEC-1  watershed 
response.  Therefore,  numerical  approxima¬ 
tions  of  the  derivatives  are  used. 


TABLE  2 


Rate  Parameter  Estimates 

Parameter  Parameter 

Initial 

Number 

Number 

Value 

1 

TC+R  /(TAREA) /TRHR 

2 

R/(TC*R) 

0.5 

3 

COEF 

0.07 

4 

STRKR 

0.2 

5 

STRKS 

0.2 

6 

RTI0K 

2.0 

7 

ERAIN 

0.5 

8 

FRZTP 

0.0 

9 

DLTKR 

0.5 

10 

RTIOL 

2.0 

TAREA  * 

Drainage  area,  in  square  miles 

TRHR  » 

Computation  interval. 

in  hours 

TABLE  3  Program  HEC-1  Default 

Initial  Routing  Parameter  Estimates 

Parameter  Initial 

Name  Value 


NSTPS 

1 

NSTDL 

1 

LAG 

1 

AMSKK 

TRHR 

X 

.2 

TSK 

3* (TRHR) 

APPLICATION  OF  THE  CALIBRATION 
CAPABILITY 

Due  to  the  varying  quantity  and  form  of  data 
available  for  precipitation-runoff  analysis, 
the  exact  sequence  of  steps  In  application 
of  the  automatic  calibration  capability  of 
HEC-1  varies  from  study  to  study.  An  often- 
used  strategy  employs  the  following  steps 
when  using  the  complete  exponential  loss 


rate  equation: 

1.  Determine  for  each  storm  selected  for  use 
in  calibration  the  base  flow  and  reces¬ 
sion  parameters  that  are  event  dependent 
and  are  not  included  in  the  set  of  param¬ 
eters  that  can  be  estimated  automatically. 
These  parameters  are  the  recession  flow 
for  antecedent  runoff  (STRTQ),  the  dis¬ 
charge  at  which  recession  flow  begins 
(QRCSN),  and  the  recession  coefficient 
that  is  the  ratio  of  flow  at  some  time 

to  the  flow  ten  time  periods  later 
(RTIOR).  These  parameters  are  illus¬ 
trated  in  Fig.  2.  The  HEC-1  Users  Manual 
(HEC,  1973)  suggests  techniques  for 
estimating  these  parameters. 

2.  For  each  storm  at  each  gage,  determine 
the  optimal  estimates  of  all  unknown  unit 
hydrograph  and  loss  rate  parameters  using 
the  automatic  calibration  feature  of 
HEC-1. 

3.  If  ERAIN  is  to  be  estimated,  select  a 
regional  value  of  ERAIN,  based  on  anal¬ 
ysis  of  the  results  of  Step  2  for  all 
storms  for  the  representative  gages. 

4.  Using  the  optimization  scheme,  estimate 
the  unknown  parameters  with  ERAIN  now 
fixed  at  the  selected  value.  Select  an 
appropriate  regional  value  of  RTIOL  if 
RTIOL  is  unknown.  If  the  temporal  and 
spatial  distribution  of  precipitation  is 
not  well  defined,  an  initial  loss, 
followed  by  a  uniform  loss  rate  may  be 
appropriate.  In  this  case,  ERAIN  *  0 
and  RTIOL  =1.  If  these  values  are  used, 
as  they  often  are  in  studies  accomplished 
at  HEC,  Steps  2,  3,  and  4  are  omitted. 

5.  With  ERAIN  and  RTIOL  fixed,  estimate  the 
remaining  unknown  parameters  using  the 
optimization  scheme.  Select  a  value  of 
STRKR  for  each  storm  being  used  for 
calibration.  If  parameter  values  for 
adjacent  basins  have  been  determined, 
check  the  selected  value  for  regional 
consistency. 

6.  With  ERAIN,  RTIOL,  and  STRKR  fixed,  use 
the  parameter  estimation  algorithm  to 
compute  all  remaining  unknown  parameters. 
DLTKR  can  be  generalized  and  fixed  if 
desired  at  this  point,  although  this 
parameter  is  considered  to  be  relatively 
event-dependent. 

7.  Using  the  calibration  capability  of  HEC-1, 
determine  values  of  TR+R  and  R/(TC+R). 
Select  appropriate  values  of  TC+R  for  each 
gage.  In  order  to  determine  TC  and  R,  an 
average  value  of  R/(TC+R)  is  typically 
selected  for  the  region. 

8.  Once  all  parameters  have  been  selected, 
the  values  should  be  verified  by  sim¬ 
ulating  the  response  of  the  gaged  basins 
to  other  events  for  which  precipitation 
and  runoff  records  are  available. 
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FORECASTING  APPLICATION 


APPLICATIONS  OF  HEC-1  AT 
THE  HEC 

The  HEC-1  optimization  scheme  has  been  used 
by  the  HEC  In  numerous  studies  for  nearly 
10  years.  The  applications  have  focused  on 
developing  frequency  curves  for  ungaged 
locations  and  on  modelling  the  Impact  of 
basin  modifications,  of  channel  Improvements, 
or  of  additional  control  measures  at  selected 
locations. 

Many  of  the  recent  applications  of  HEC-1 
accomplished  at  the  Hydrologic  Engineering 
Center  have  employed  the  automatic  calibra¬ 
tion  scheme  In  development  of  data  for  ungaged 
areas.  Typically  in  these  studies  data  are 
available  from  stream  and  precipitation 
measurement  stations  in  the  proximity  of  a 
location  for  which  detailed  stage  and  dis¬ 
charge  data  are  unavailable  but  are  desired. 

The  automatic  calibration  technique  is  used 
to  estimate  unit  hydrograph,  loss  rate,  and 
routing  parameters  for  the  gaged  locations, 
and  this  data  is  “transferred"  to  the  ungaged 
locations  using  regression  techniques.  The 
particular  strategy  for  estimating  parameters 
and  the  methods  for  transferring  the  param¬ 
eters  to  ungaged  locations  is  a  function  of 
the  basin  characteristics,  the  available  data, 
the  parameters  that  are  found  via  calibration, 
and  the  time  and  money  available  for  the  study. 

While  the  sequence  of  steps  for  estimation  of 
all  parameters  of  this  ralnfal 1 -runoff  model 
has  been  employed  in  at  least  one  major  study, 
the  flexibility  gained  by  use  of  four  param¬ 
eters  in  the  exponential  loss  rate  equation 
Is  not  always  necessary.  Often  ERAIN  is  set 
to  zero,  RTIOL  is  set  to  one,  and  calibration 
proceeds  with  Step  5  to  this  sequence.  This 
approach  has  been  employed  in  studies  of  the 
Shell  pot  and  Naaman  Creeks  (HEC,  1976),  the 
Schuylkill  River  (HEC,  1976),  the  Maurice 
River  (HEC,  1976),  and  the  Lehigh  River 
(HEC,  1978).  In  a  hydrologic  study  on  the 
Oconee  River  Basin,  ERAIN  and  DLTKR  were  set 
to  zero,  yielding  a  simple  exponential  decay 
loss  function  (1976). 

The  Soil  Conservation  Service's  dimensionless 
unit  graph  and  rainfall -runoff  relationship 
based  on  the  soil  classification  curve  numbers 
were  recently  added  to  HEC-1.  An  application 
was  made  on  the  Pennypack  Creek  as  part  of  an 
expanded  flood  plain  Information  report 
(HEC,  1978).  The  145  km*  study  area  was 
broken  Into  65  subareas,  each  requiring  unit 
graph,  loss  rate,  and  channel  storage  param¬ 
eters.  The  curve  number  (CN)  and  lag  param¬ 
eters,  which  are  the  only  variables  necessary 
to  define  the  SCS  unit  graph  and  loss  rate, 
were  estimated  via  optimization  for  the  gaged 
basins.  These  were  used  as  a  guide  in  estab¬ 
lishing  CN  and  LAG  values  for  the  ungaged  sub- 
areas  . 


In  addition  to  the  previously  described 
applications  of  HEC-1,  the  model  has 
recently  been  applied  to  develop  reservoir 
inflow  forecasts  for  W.  Kerr  Scott  Reservoir 
on  the  Yadkin  River  of  North  Carolina.  For 
this  application,  the  basic  model  was 
modified  so  that  the  calibration  technique 
could  be  used  to  update  sequentially  the 
model  parameters.  The  parameters  are  then 
used  to  calculate  forecasted  streamflows. 

In  the  W.  Kerr  Scott  system,  20  gages,  nine 
of  them  recording,  were  available  to  deter- 2 
mine  mean  areal  precipitation  in  the 900  knr 
basin.  Seven  storm  events  were  modelled, 
using  the  optimization  procedure  to  estimate 
TC,  R,  ERAIN,  RTIOL,  and  STRKR  for  the  basin. 
Conceptual ly,  these  parameters  will  remain 
constant  from  storm  to  storm  on  the  same 
basin,  with  DLTKR  alone  indicating  the  ante¬ 
cedent  basin  wetness.  In  reality,  all  the 
parameters  vary  due  to  storm  centering,  to 
inaccuracies  in  precipitation  data,  to  non¬ 
homogeneity  of  the  basin,  and  to  the  approx¬ 
imate  nature  of  the  hydrologic  model.  Three 
combinations  of  estimated  and  fixed  parameters 
were  investigated.  In  the  first  case,  basin 
average  values  of  TC  and  R  were  used,  and 
the  four  loss  rate  parameters  were  established 
by  calibration.  In  the  second  case,  ERAIN  was 
set  to  zero  and  RTIOL  was  set  to  one,  while 
the  remaining  parameters  were  estimated  by 
calibration.  Finally,  all  six  parameters 
were  estimated  via  calibration.  When  tested 
for  seven  storm  events,  the  latter  approach 
provided  more  accurate  predictions  when  the 
time  of  forecast  occurred  before  the  peak 
of  the  hydrograph.  In  each  case,  the  param¬ 
eters  were  estimated  using  the  computed  re¬ 
servoir  inflows  and  observed  mean  areal 
precipitation  that  were  available  at  the  time 
of  forecast.  The  purpose  of  the  calibration 
was  to  adjust  loss  rates  to  reflect  the  ante¬ 
cedent  conditions  and  fine  tune  the  unit  graph 
parameters  to  best  fit  each  storm. 

Forecast  Procedure 

In  order  to  test  the  procedure  in  a  realistic 
forecasting  application,  data  from  four  storms 
were  used  in  a  situation  that  is  similar  to 
that  encountered  by  Corps'  field  offices. 
Whenever  the  calculated  reservoir  inflow 
exceeded  50  cms  (twice  the  long  term  average 
inflow  for  Scott  Reservoir),  a  forecast  was 
calculated  using  the  data  available  at  that 
time.  Subsequent  forecasts  were  calculated 
every  two  hours,  using  additional  data  obser¬ 
ved  between  forecasts,  until  the  Inflow 
hydrograph  started  to  recede.  In  addition, 
forecasts  were  Issued  whenever  subsequent 
rises  occurred  during  the  receding  limb. 

Thus  forecasts  were  made  for  all  the  peaks 
in  a  complex  hydrograph. 
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Two  strategies  were  investigated  for  estab¬ 
lishing  the  initial  values  for  calibration 
of  the  unit  graph  and  loss  rate  parameters. 


The  initial  values  for  TC,  R,  ERAIN,  RTIOL 
for  method  1  were  based  on  average  values 
determined  from  the  calibration  of  seven 
storms  while  program  default  values  were  used 
for  STRKR  and  DLTKR.  These  starting  values 
were  used  for  every  forecast  during  the  event, 
and  did  not  consider  the  results  of  the  pre¬ 
vious  forecast.  Method  2  employed  the  same 
starting  values  for  the  first  forecast  of 
each  storm,  but  initial  values  for  subsequent 
forecasts  were  set  equal  to  optimal  parameter 
values  from  the  previous  forecast. 

Analysis  of  Forecasting  Results 

The  performance  of  a  flood  forecast  model  can 
be  fully  evaluated  only  after  the  event  has 
occurred  and  the  complete  observed  inflow 
hydrograph  is  available.  For  reservoir 
operation,  both  the  shape  and  the  volume  of 
the  flood  hydrograph  are  important.  The 
following  statistical  measures  were  defined 
to  measure  discrepancies  between  the  fore¬ 
casted  and  observed  inflows,  beginning  when 
the  forecast  is  issued  and  ending  when  the 
observed  flow  recedes  to  20%  of  the  peak  flow: 


Volume  Error  - 

N  N 

^  (QS,  -  QO.)/  2  QO. 
i*l  1  1  1-1  1  (9) 

Average  Discharge  Error  = 

N 

1/N  s  (IQS.  -  QOJ/QOJ  *  100 

1»1  1  11  (10) 


where  QO.  3  observed  discharge  at  ordinate  i; 
QS^  =*  simulated  discharge  at  ordinate  i;  N  - 
number  of  hydrograph  ordinates. 

The  overall  accurancy  of  the  forecasts  is'  a 
function  of  storm  type  and  intensity,  repre- 
sentati  veness  of  the  computed  basin  average 
precipitation,  and  the  degree  of  storm  develop¬ 
ment,  Despite  the  difficulties  in  providing 
high  quality  data  for  the  calibration,  the 
forecasted  volumes  for  the  four  storms  are 
reasonable;  the  predicted  inflow  volumes  are 
within  25  percent  of  the  observed  values  8 
to  10  hours  prior  to  the  peak.  As  expected, 
forecasts  made  before  the  total  precipitation 
had  fallen  generally  underestimated  the  total 
runoff  volume. 

In  order  to  Improve  the  performance  of  the 
forecasting  scheme,  the  sources  of  differences 
between  the  observed  and  forecasted  discharge 
must  be  Isolated,  Typically  these  sources 
Include  errors  inherent  within  the  forecast¬ 
ing  procedure,  errors  in  the  basin  hydrologic 
model,  and  errors  in  the  model  calibration. 

Errors  inherent  within  any  forecasting  pro¬ 
cedure  include  uncertainty  in  measurement  of 
precipitation  and  streamflow  data,  and  under¬ 


estimation  of  the  hydrograph  before  the  total 
precipitation  occurs.  In  the  Yadkin  applica¬ 
tion,  inflow  into  the  W.  Kerr  Scott  Reser¬ 
voir  was  calculated  using  the  change  in  eleva¬ 
tion  of  the  pool  level,  elevation-storage 
tables,  and  the  reservoir  releases  determined 
by  the  gate-discharge  rating  tables.  Fluctu¬ 
ations  in  the  reservoir  pool  due  to  wind  set¬ 
up  and  seiching,  combined  with  an  unknown 
amount  of  gate  slippage  caused  the  calculated 
inflows  to  oscillate  considerably .  The  adopted 
inflow  hydrograph  used  for  calibrating  and 
forecasting  was  based  on  a  smoothed  curve,  so 
considerable  errors  could  occur  in  the  slope 
of  the  rising  limb.  Because  this  early  por¬ 
tion  of  the  observed  hydrograph  is  used  to 
verify  and  adjust  the  simulated  hydrograph, 
the  calibration  scheme  could  have  forced  a 
fit  with  erroneous  data  and  thus  could  have 
degraded  the  forecasted  hydrograph. 

The  basic  structure  of  the  hydrologic  model 
and  the  manner  it  is  applied  can  also  lead  to 
errors.  Because  of  the  appreciable  lag  time 
between  the  measurement  of  rainfall  on  the 
basin  and  its  occurrence  at  the  streamgage, 
when  a  forecast  is  issued  early  in  the  flood 
event,  only  a  small  portion  of  the  recorded 
precipitation  may  contribute  to  the  simulated 
discharge  prior  to  the  time  of  forecast. 

When  this  occurs  it  is  possible  that  the  optim¬ 
ization  scheme  may  cause  substantial  changes 
in  the  model's  parameters  without  a  corre- 
esponding  change  in  the  fit  between  the  sim¬ 
ulated  and  observed  discharge.  This  effect 
was  evident  when  alternative  initial  param¬ 
eter  values  were  employed.  In  17  of  25 
forecasts,  method  2  exhibited  a  lower  cal¬ 
ibration  error  (error  prior  to  forecast), 
while  in  all  but  two  cases  method  1  exhibited 
a  lower  total  error  for  the  event.  In  seven 
of  the  cases  for  which  method  2  had  the  lower 
calibration  error,  the  total  event  error  was 
two  to  five  times  greater  than  that  with 
method  1.  This  demonstrates  that  an  excel¬ 
lent  fit  between  the  observed  and  simulated 
discharges  does  not  necessarily  lead  to  a 
correspondingly  good  forecast. 

Several  characteristics  of  the  parameter 
estimation  scheme  can  influence  the  amount 
of  forecast  error.  First,  if  the  unit  hydro¬ 
graph  and  loss  rate  parameters  are  constrained 
too  loosely,  as  in  the  HEC-1  scheme, 
the  adopted  parameters  may  be  unreasonable. 
Method  2,  which  consistently  yielded  better 
results  during  calibration,  often  calculated 
unrealistic  values  of  TR,  R  and  RTIOL.  A 
second  factor  is  the  allowable  error  between 
the  observed  and  simulated  values  in  the 
calibration  step.  To  reduce  the  calibration 
error,  the  parameter  estimation  algorithm 
may  produce  significant  fluctuations  in  par¬ 
ameters  when  only  a  small  portion  of  precip¬ 
itation  and  discharge  data  are  available. 
However,  if  the  allowable  tolerance  is  in¬ 
creased  sufficiently  early  in  the  storm  event 
and  is  gradually  reduced,  the  number  and  mag¬ 
nitude  of  parameter  modifications  would  be 
reduced.  This  concept  was  adopted  for  a  flood 
forecasting  model  developed  by  the  National 
Weather  Service  (1979). 
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In  suanary,  the  performance  of  the  calibra¬ 
tion  scheme  was  only  one  contributing  factor 
to  the  accuracy  of  the  flood  forecast.  While 
the  addition  of  a  variable  tolerance  or  more 
limiting  constraints  to  the  present  method 
may  be  beneficial,  the  results  of  the  HEC-l 
model  were  quite  satisfactory. 


ALTERNATIVE  OPTIMIZATION  ALGORITHMS 
FOR  CALIBRATION 

In  1978  the  Hydrologic  Engineering  Center  began 
a  study  of  alternative  techniques  for  estima¬ 
tion  of  parameters  for  HEC-l.  This  study  was 
motivated  by  increasing  use  by  Corps  field 
offices  of  the  automatic  calibration  feature 
of  the  model.  In  these  many  applications  of 
the  model,  any  significant  improvement  in  the 
parameter  estimation  technique  is  desirable, 
in  terms  of  a  reduction  in  the  cost  of  cal¬ 
ibration  or  in  terms  of  a  reduction  in  the 
error  of  the  parameter  estimates. 

Analysis  of  the  existing  univariate  search 
technique  employed  for  parameter  estimation 
in  HEC-l  indicated  that  the  technique  often 
does  not  effectively  handle  constraints  on 
the  parameters,  so  a  nonlinear  programming 
algorithm  that  does  was  sought.  An  addi¬ 
tional  restriction  on  selection  of  an  op¬ 
timization  scheme  is  imposed  by  the  number 
of  parameters  that  must  be  estimated.  If  all 
parameters  of  the  precipitation-runoff  process 
are  to  be  estimated,  the  optimization  algo¬ 
rithm  must  be  capable  of  solving  efficiently 
a  problem  of  as  many  as  ten  variables. 

The  alternative  optimization  technique  sel¬ 
ected  for  initial  trial  application  was  a 
version  of  the  random  search  method  of  Box 
(1965),  as  suggested  by  Johnston  and  Pilgrim 
(1976),  Chu  and  Bowers  (1978),  and 
Sorooshian  and  Dracup  (1978).  This  technique 
is  an  extension  of  the  polyhedron  search  of 
Nelder  and  Mead  (1965),  and  as  programmed  at 
HEC,  uses  the  precipitation-runoff  model 
directly  in  evaluation  of  the  objective  func¬ 
tion.  Implementation  of  the  technique  re¬ 
quires  definition  of  explicit  upper  and  lower 
bounds  on  the  variables;  this  proved  to  be  a 
formidable  task  because  the  parameters  of  the 
rafnfal 1 -runoff  model  are  not  related  clearly 
to  physical  attributes  of  the  drainage  basin. 
Currently  the  constraints  are  defined  on  the 
basis  of  knowledge  of  the  limitations  of  the 
mathematical  functions  combined  with  exper¬ 
ience  in  manual  calibration  of  the  model  for 
Corps'  studies  nation-wide. 

Although  the  investigation  of  alternative 
optimization  techniques  or  calibration  schemes 
is  not  yet  complete,  preliminary  analyses 
indicate  the  following:  (1)  the  Box  tech¬ 
nique  does  not  select  an  optimal  set  of  param¬ 
eter  values  in  significantly  less  time  or  with 
fewer  function  evaluations  than  does  the 
univariate  gradient  algorithm  currently 
employed;  and  (2)  the  Box  technique  does  not 
select  parameter  values  that  yield  a  sig¬ 
nificantly  better  reconstitution  of  runoff 


events,  as  measured  by  the  least-squares 
objective  function.  As  a  result  of  the 
research,  the  following  conclusions  were 
reached  regarding  the  parameter  estimation 
procedure:  (1)  for  efficiency,  an  automatic 
calibration  technique  should  allow  specifi¬ 
cation  of  and  should  exploit  the  availability 
of  appropriate  initial  values  of  the  param¬ 
eters;  and  (2)  the  least-squares  objective 
function  is  relatively  insensitive  to  vari¬ 
ations  of  the  parameters  of  HEC-l. 

The  Box  technique  is  executed  by  initially 
evaluating  the  objective  function  for  sets 
of  parameter  values  scattered  randomly 
throughout  the  feasible  region  (a  "complex"). 
Box  recommends  that  2n  sets  of  values  be 
established  for  the  initial  complex  where  n 
equals  the  number  of  unknown  parameters.  For 
ten  parameters,  this  requires  20  evaluations 
of  the  objective  function,  each  requiring 
simulation  of  watershed  response.  The  uni¬ 
variate  technique,  with  certain  heuristic 
rules  for  dealing  with  violated  constraints 
and  with  good  initial  estimates  often  selects 
near-optimal ,  acceptable  parameter  estimates 
with  approximately  the  same  number  of  func¬ 
tion  evaluations.  Furthermore,  the  random 
search  techniques  does  not  exploit  the 
knowledge  of  good  initial  estimates  of  the 
parameters.  Experience  at  HEC  indicates  that 
this  is  critical  because  these  estimates  often 
reduce  substantially  the  effort  to  calibrate 
and  because  automatic  calibration  schemes  may 
chose  an  unreasonable,  false  optimum  other¬ 
wise. 

The  least-squares  objective  function,  although 
widely  accepted  for  application  in  model 
calibration,  in  many  cases  is  insensitive  to 
variations  in  the  parameters  of  HEC-l  and 
causes  premature  termination  of  the  search. 
This  difficulty  is  related  also  to  interaction 
among  the  variables  of  the  watershed  model. 
However,  acceptance  of  any  modifications  to 
the  model  is  not  likely,  so  alternative  objec¬ 
tive  functions  are  being  evaluated.  These 
functions  include  those  suggested  by  Manley 
(1978)  and  other  functions  suggested  through 
research  at  HEC. 

Additional  tasks  to  be  completed  in  the  area 
of  parameter  estimation  for  the  HEC-l  water¬ 
shed  model  include  the  following:  (1)  pro¬ 
gramming  and  testing  other  nonlinear  program¬ 
ming  algorithms  for  parameter  estimation; 

(2)  programming  and  testing  alternative 
objective  functions  for  calibration;  (3) 
further  comparing  the  existing  parameter 
estimation  algorithm  with  the  alternative 
algorithms;  and  (4)  applying  the  knowledge 
gained  from  parameter  estimation  research 
with  program  HEC-l  to  other  programs  devel¬ 
oped  and  supported  by  the  HEC. 


SUMMARY 

To  satisfy  the  need  for  a  precipitation- 
runoff  model  for  application  in  water  resour¬ 
ces  planning  and  management  by  the  Corps  of 
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Engineers,  computer  program  HEC-1  was  devel¬ 
oped.  This  model  Includes  algorithms  to 
accomplish  the  following  tasks  necessary  to 
simulate  watershed  response: 

1.  Determine  effective  precipitation. 

2.  Compute  the  subarea  runoff  due  to  the 
effective  precipitation. 

3.  Route  and  combine  the  subarea  runoff 
hydrographs , 

In  addition,  HEC-1  includes  the  capability  to 
determine  automatically  the  parameters  of  the 
functions  employed  in  the  simulation.  This 
is  accomplished  using  Newton's  technique  to 
minimize  a  weiqhted  least- squares  objective 
function.  Currently,  alternative  optimi¬ 
zation  techniques  and  alternative  objective 
functions  are  being  evaluated. 

The  parameter  estimation  capability  of  pro¬ 
gram  HEC-1  has  been  employed  in  a  variety 
of  studies  at  the  Hydrologic  Engineering 
Center.  These  applications  have  focused  on 
modelling  the  impact  of  basin  modifications, 
of  channel  Improvements,  of  various  flood- 
control  measures,  and  on  developing  fre¬ 
quency  curves  for  ungaged  watersheds. 

The  parameter  estimation  technique  of  HEC-1 
has  been  extended  recently  to  update  sequen¬ 
tially  parameter  estimates  for  flood  fore¬ 
casting.  In  these  applications,  computed 
reservoir  inflows  and  observed  mean  areal 
precipitation  available  at  the  time  of  fore¬ 
cast  are  used  to  estimate  model  parameters . 
These  parameters  are  used  then  to  estimate 
future  reservoir  inflows,  using  the  simula¬ 
tion  capability  of  the  program.  The  results 
of  the  applications  are  satisfactory  for 
application  to  flood-control  reservoir  oper¬ 
ation. 
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